Accurate calculation of the local density of optical states in 
inverse-opal photonic crystals 



Ivan S. Nikolaev'^, Willem L. Vos^'^, and A. Femius Koenderink^ * 

^ Center for Nanophotonics, FOM Institute for Atomic and Molecular Physics (AMOLF), Kruislaan 407, NL-1098SJ 

Amsterdam, The Netherlands 



o 
o 

(N 
D 



C/3 

o 

> 

q : 
o ■ 

• I— I ■ 

^ : 

Oh. 



> 

o 
in 

00 

O 
On 
O 



X 



'Complex Photonic Systems, MESA'^ Institute for Nanotechnology, University of Twente, The Netherlands 

Corresponding author: f.koenderink@amolf.nl 

Submitted to OS A Doc 23, 2008 

We have investigated the local density of optical states (LDOS) in titania and silicon inverse opals — 
three-dimensional photonic crystals that have been realized experimentally. We used the H-field plane-wave 
expansion method to calculate the density of states and the projected local optical density of states, which 
are directly relevant for spontaneous emission dynamics and strong coupling. We present the first quantitative 
analysis of the frequency resolution and of the accuracy of the calculated local density of states. We have 
calculated the projected LDOS for many different emitter positions in inverse opals in order to supply a 
theoretical interpretation for recent emission experiments and as reference results for future experiments and 
theory by other workers. The results show that the LDOS in inverse opals strongly depends on the crystal 
lattice parameter as well as on the position and orientation of emitting dipoles. 
OCIS codes: 160.5298,260.2510,270.5580,290.4210 



1. Introduction 

Photonic crystals are metamaterials with periodic varia- 
tions of the dielectric function on length scales compara- 
ble to the wavelength of light. These dielectric compos- 
ites are of a keen interest for scientists and engineers 
because they offer exciting ways to manipulate pho- 
tons [1,2]. Of particular interest are three-dimensional 
(3D) photonic crystals possessing a photonic bandgap, 
i.e. a frequency range where no photon modes exist at 
all. Photonic bandgap materials possess a great potential 
for drastically changing the rate of spontaneous emission 
and for achieving localization of light [3-6] . Control over 
spontaneous emission is important for many applications 
such as miniature lasers [7] , light-emitting diodes [8] and 
solar cells [9]. According to Fermi's golden rule, the rate 
of spontaneous emission from quantum emitters such as 
atoms, molecules or quantum dots is proportional to the 
'local radiative density of states' (LDOS) [10,11], which 
counts the number of electromagnetic states at given fre- 
quency, location and orientation of the dipolar emitters. 
In addition, nonclassical effects can occur that are be- 
yond Fermi's Golden Rule [2,11]. These strong coupling 
phenomena, such as fractional decay, rely on coherent 
coupling of the quantum emitter to a sharp feature in 
a highly structured LDOS. Whether these effects are in 
fact observable in real photonic crystals depends on how 
rapid the LDOS changes within a small frequency win- 
dow [12]. Therefore local density of states calculations 
for experimentally realized photonic crystals are essen- 
tial to assess current spontaneous experiments and fu- 
ture strong coupling experiments alike, provided that the 
calculations are accurate and have a controlled frequency 
resolution. 



LDOS effects on spontaneous emission in photonic 
crystals have been experimentally demonstrated in a 
variety of systems. Since only three-dimensional (3D) 
crystals promise full control over all optical modes with 
which elementary emitters interact, many groups have 
pursued 3D photonic crystals. Fabrication of such pe- 
riodic structures with high photonic strengths is, how- 
ever, a great challenge [2,13]. Inverse opals are among 
the most strongly photonic 3D crystals that can be fabri- 
cated relatively easily using self-assembly methods. Such 
crystals consist of fee lattices of close-packed air spheres 
in a backbone material with a high dielectric constant 
[14-18]. In these inverse opals, continuous- wave exper- 
iments on light sources with a low quantum efficiency 
revealed inhibited radiative emission rates [19,20]. En- 
hanced and inhibited time-resolved emission rates have 
been observed from highly-efficient emitters in 3D in- 
verse opals [21]. Simultaneously, several groups have re- 
alized that emission enhancement and partial inhibition 
can also be obtained in 2D slab structures [22-26] . 

Following the first calculations by Suzuki and Yu [27], 
several papers have reported on calculations of the lo- 
cal density of optical states in photonic crystals using 
both time domain [28-30] and frequency domain meth- 
ods [31-34]. Unfortunately, analysis of experimental data 
for emitters in 3D crystals in terms of these LDOS calcu- 
lations is compounded by several problems in literature: 

1. Most prior calculations were performed for model 
systems that don't correspond to structures used in 
experiments, and for emitter positions that are not 
probed in experiments. 

2. The accuracy of the reported LDOS has never been 
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discussed, hampering comparison to experiments. 

3. The frequency resolution of the reported LDOS has 
remained unspecified. Therefore sharp features of 
relevance for non-classical emission can not be as- 
sessed [11,12]. 

4. Many previously reported LDOS calculations are er- 
roneous for symmetry reasons, as pointed out by 
Wang et al. [33]. 

In this paper we aim to overcome all these problems. 
We benchmark the accuracy and frequency resolution for 
our LDOS calculations to allow quantitative comparison 
with experiments and with quantum optics strong cou- 
pling requirements. Since prior LDOS calculations are 
scarce (and partly erroneous, item 1 above) we present 
sets of LDOS's for experimentally relevant structures, 
and for spatial positions where sources can be practi- 
cally placed. Specifically, we model the spatial distribu- 
tion of the dielectric function e(r) in such a way that it 
closely resembles e(r) in titania (Ti02) [14, 15] and sili- 
con (Si) inverse-opal photonic crystals [17,18]. For these 
two structures we calculated the LDOS at various posi- 
tions in the crystal unit cell and for specific orientations 
of the transition dipoles. The results on the Ti02 inverse 
opals are relevant for interpreting recent emission exper- 
iments [21,35]. To aid other workers to interpret their ex- 
periments and to benchmark their codes, we make data 
sets that we report in figures throughout this manuscript 
available as online material. 

The paper is arranged as follows: in Section [2l we 
present a detailed description of the method by which 
we have calculated the photonic band structures and the 
LDOS. We discuss the accuracy and frequency resolution 
of our calculations. In Section [3] we compare our compu- 
tations with the known DOS in vacuum and with earlier 
results on the DOS and LDOS [33] in 3D periodic struc- 
tures. SectionlUdescribes the LDOS in inverse opals from 
Ti02, and in Section [5] we present results of the LDOS 
in Si inverse opal photonic band gap crystals. 

2. Calculation of local density of states 

A. Introduction 

The local radiative density of optical states is defined as 

7V(r,a;,ed) = — i-r V / dk(5(w- cj„,k)|ed •E„,k(r)l^, 

(1) 

where integration over k-vcctor is performed over the 
first Brillouin zone, n is the band index and e^ is the 
orientation of the emitting dipole. The total density of 
states (DOS) is the unit-cell and dipole-orientation av- 
erage of the LDOS defined as N(u:) = J^^ dkS{uj — 
'j-'n,k)- The important quantities that determine the 
LDOS are the eigenfrequencies LUn^k and electric field 
eigenmodes E„.k(r) for each k- vector. Calculation of 
these parameters will be discussed below in Section |B] 



The expression for the LDOS contains the term je^; • 
En,k(r)l that depends on the dipole orientation e^. It is 
important to realize that in photonic crystals the vec- 
tor fields E„_k(r) are not invariant under the lattice 
point-group operations a, as first reported in Ref. [33]. 
Explicitly, this means that the projection of the field 
of a mode at wave vector k on the dipole orientation 
e^ (i.e. \ed ■ E„.k(r)l)is not identical to the projection 
je^ • E„ (j[k] (r) 1 of the symmetry related modes with wave 
vectors a[k] on the same e^. As a consequence one can 
not calculate the LDOS for a specific dipole orientation 
by restricting the integral in Eq. ^ over the irreducible 
part (l/48th) of the Brillouin zone, since symmetry re- 
lated wave vectors do not give identical contributions. 
Unfortunately, in many previous reports on the LDOS, 
this reduced symmetry for vector modes as compared to 
scalar quantities was overlooked, resulting in erroneous 
results [33]. In general, the only symmetry that can be 
invoked to avoid using the full Brillouin zone for LDOS 
calculations is inversion symmetry, which corresponds to 
time reversal symmetry. Consequently, correct results re- 
quire that exactly half of the Brillouin zone is considered 
for LDOS calculations, rather than the irreducible part 
of the Brillouin zone that was used in most previous lit- 
erature. We have explicitly verified that our implemen- 
tation (using the > half of the Brillouin zone) re- 
sults in the same LDOS on symmetry related positions, 
provided that one also takes the concomitant symmetry 
related dipole orientation into account. Furthermore, our 
calculations confirm the claim by Wang et. al. that this 
required symmetry is only recovered upon integration 
over half the Brillouin zone, rather than over just the 
irreducible part as considered in, e.g., Ref. [31,32]. 

B. Plane-wave expansion 

We use the H-ficld inverted plane wave expansion 
method [31,36,37] to solve for the electromagnetic field 
modes in photonic crystals. For nonmagnetic materials, 
it is most convenient to solve the wave equation for the 
H(r) field [38] 

V X [e(r)-iV X H(r)] = ^^H(r). (2) 

because the operator V x e(r)^^Vx is Hermitian, and 
consequently has real eigenvalues w^/c^ [1,31,36,37]. Be- 
cause of the periodicity of the dielectric function e(r) in 
photonic crystals, the field modes Hk(r) of the eigen- 
value problem Eq. ([2]) satisfy the Bloch theorem [39] : 

Hk(r) = e''''-Uk(r). (3) 

These Bloch modes arc fully described by the wavevector 
k and the periodic function Uk(r), which has the period- 
icity of the crystal lattice so that Uk(r) = Uk(r -I- R). To 
solve the wave equation, the inverse dielectric function 
and the Bloch modes are expanded in a Fourier series 
over the reciprocal-lattice vectors G: 

e(r)-i ^ rj{r) = J2 me'^'" and (4) 

G 
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k „i(k+G)-r 



Une 



(5) 



where tjg and Uq are the 3D Fourier expansion coefS- 
cients of respectively 77(r) and Uk(r). Substituting these 
expressions into the H-field wave equation in Eq. ([2]), we 
obtain a hnear set of eigenvalue equations: 



J2 m-G'ik + G) X [(k + G') X u^!^] 

G' 



^nW „»,k 
„2 G • 



(6) 

This infinite equation set with the known parameters G 
and riG-G' determines all allowed frequencies w„(k) for 
each value of the wave vector k, subject to the transvcr- 
sality requirement V • Hk(r) — 0. Due to the periodicity 
of Uk(r), we can restrict k to the first Brillouin zone. For 
each wave vector k, there is a countably infinite num- 
ber of modes with discretely spaced frequencies. All the 
modes are labeled with the band number n in order of 
increasing frequency and are described family of 

continuous functions a;„ (k) of k. 

To compute the eigenfrequencies a;„(k) and the ex- 
pansion coefficients of the eigenmodes Uq'', the infinite 
equation set is truncated. By restricting the number of 
reciprocal-lattice vectors G to a finite set G with Nq ele- 
ments, Eq. ^ is limited to a SNg dimensional equation 
set. In our implementation we choose the truncated set G 
to correspond to the set of all reciprocal lattice vectors 
within a sphere centered around the origin of k-spacc. 
The transversality of the H-ficld gives an additional con- 
dition on the eigenmodes: (k + G) •u^'' = 0, which elim- 
inates one vector component of u^*'. Following Ref. [31], 
for each k -f G one needs to find two orthogonal unit vec- 
tors eJj.'^Q that form an orthogonal triad with k + G. By 
expressing the cigenmode expansion coefficients in the 
plane normal to k + G as Uq'' = "c'l^k+G +^G'2^k+Gi 
we remove one third of the unknowns. Then, Eq. ([6]) be- 
comes 



J2 »7G-G'|k + G||k- 
G'ee 



G'l- 



^k+G ^k+G' 
'^k+G *^k+G' 


•^k-l-G 
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n,k 
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VGeg. 



(7) 

To find the matrix of Fourier coefficients tjg-G' , we used 
the method of Refs. [31,36]. The coefficients are com- 
puted by first Fourier-transforming the dielectric func- 
tion e(r), and then truncating and inverting the resulting 
matrix. As first noted by Ho, Chan and Soukoulis [36], 
using the inverse of eg-G' rather than 7?g-G' dramati- 
cally improves the (poor) convergence of the plane- wave 
method that is associated with the discontinuous nature 
of the dielectric function [37]. A rigorous explanation 
for this improvement was put forward by Li [40], who 



studied the presence of Gibbs oscillations in the trun- 
cated Fourier expansion of products of functions with 
complementary jump discontinuities. Using the H-field 
inverted matrix plane wave method, the frequencies ob- 
tained with Ng = 725 (for fee structures) deviate by 
less than 0.5 % from the converged band structures [31]. 
Solving Eqs. ([7]) gives the frequencies ti;„(k) and the 
Fourier expansion coefficients for the H-field eigenmodes 
H„^k(r) needed to calculate the LDOS in the photonic 
crystal. The required E-ficlds E„,k(r) are obtained using 
the Maxwell equation dG/dt = V x H: 



E„.k(r) 



1 



^«(k)eo 



G,G'ee 



n.k 2 
"G.l^k+G 



ri.k 1 
^G,26k+G 



,l(k+G')T 



(8) 

From the orthonormality of the eigenvectors of Eq. ((T]) 
it follows that the Bloch functions H„.k(r) and E„^k(r) 
defined above satisfy the orthonormality relations: 



H„,k(r)-H;, k'(r)rfr-<5(k-k')<5„ 



(9) 



BZ 



I e(r)E„,k(r) • E;,.k, (r)dr = ,5(k - k')5„,„- . (10) 

JBZ 

It should be noted in the definition of E„.k that the mul- 
tiplication by l/e(r) to calculate E from D is not in the 
denominator in front of the Fourier expansion. Rather it 
appears as a matrix multiplying the D-field in Fourier 
space, i.e., within the sum over reciprocal lattice vec- 
tors. This ordering ensures that complementary jumps 
in D and l/e(r) cancel, even for the truncated Fourier 
series, as can be easily checked by plotting calculated 
mode profiles for TM modes in 2D crystals. 

C. Frequency resolution and aeeuraey of LDOS 

We are not aware of any previous report that bench- 
marks the accuracy of the calculated local radiative den- 
sity of states or that specifies the frequency resolution. 
Motivated by the requirements for accurate results for 
comparison with experiments and for judging the utility 
of crystals for non-classical emission dynamics, we con- 
sider the accuracy and resolution of our approximation 
to the LDOS integral in Eq. ([1]). The main approxima- 
tion is to replace the integration over wave vector k by 
an appropriately weighted summation over a discrete set 
of wave vectors on a discretization grid [41]. Either in- 
terpolation schemes [31] or simple histogramming ('root- 
sampling') methods are used to compute the LDOS. The 
k-point grid density sets the number of k-points in the 
Brillouin zone, N)^. The accuracy of the resulting LDOS 
approximation is set by the the density of grid points 
that is used to discretize the wave vector integral. Due 
to the transparent relation between the accuracy of the 
DOS, the frequency resolution and the k- vector sampling 
resolution, we will focus on the simple histogramming 
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Fig. 1. (Color online) DOS per volume in units for 
vacuum modeled as an empty fee crystal. The calculated 
DOS shown by red histogram bars is compared to the 
analytically derived behavior (curve) . In vacuum the 
DOS per volume equals the dipole-averaged LDOS. 




Aca/Ak [in units of c] 

Fig. 2. Average absolute deviation of the calculated DOS 
from the exact total DOS N{uj) for an 'empty crystal.' 
The average runs over the frequency range < a; < 
I'Kcja and the deviation is in units of the DOS N{lo) 
per volume at Lua/2Trc = 1, i.e., in units 4/a^c. In accor- 
dance with Eq. (jlip . the error is inversely proportional 
to the ratio Aw/Afc of the histogram bin width Auj to 
the integration grid spacing Afc. Symbols correspond to 
integration using the number of k-points Nk = 280, 770, 
1300, 2480, 2992 and 3570 (□, ■, o, •, 0, ♦) respectively, 
with various Atu. 



approach. For the LDOS computations one chooses a 
certain frequency bin width Alo to build an LDOS his- 
togram. For a desired frequency resolution Alu and a 
desired accuracy for the LDOS content N{r, w, Bj;) Aw in 
each frequency bin, one needs to choose an appropriate 
k-vector spacing Afc that depends on the steepness of 
the sampled dispersion relation w(k;). Indeed, the useful 
frequency resolution Alu of a histogram of the DOS and 
LDOS is limited by the resolution Afc of the grid in k 



space to be 

Aiv Ak\Vi,uj\, (11) 

as detailed in Ref. [42] for the electronic DOS. This cri- 
terion relates the separation between adjacent k-grid 
points to their approximate frequency spacing via the 
group velocity. If histogram bins are chosen too narrow 
compared to the expected frequency spacing between 
contributions to the discretized LDOS integral, unphysi- 
cal spikes appear in the approximation, especially in the 
limit of small w, where the group velocity IVkw] is usu- 
ally largest. Apart from full gaps in the LDOS, photonic 
crystals also promise sharp lines at which the LDOS is 
enhanced, which are important for non-classical emis- 
sion dynamics. Hence it is especially important to dis- 
tinguish sharp spikes that are due to histogram binning 
noise from true features. Unfortunately, many reports in 
literature feature sharp spikes that are evidently binning 
noise (wave vector undersampling errors), as they occur 
in the long wavelength limit, below any stop gap. 

To improve the resolution without adding time- 
consuming diagonalizations, several interpolation 
schemes have been suggested [42]. Within the his- 
togramming approach, an interpolation scheme essen- 
tially improves binning statistics by adding histogram 
contributions from intermediate grid points on the 
assumption that quantities vary linearly between grid 
points. While interpolation decouples the binning noise 
from the frequency bin width Alo, resulting in arbitrarily 
smooth LDOS approximations, it has the disadvantage 
of obscuring the intrinsic relation between frequency 
resolution and wave vector resolution in Eq. (jlip . 

A good benchmark for the binning noise of the k-space 
integration, independent of the convergence of the plane- 
wave method, is to calculate the LDOS or DOS of an 
'empty' crystal, with uniform dielectric constant equal 
to unity. Such an 'empty' crystal represents a limit of 
zero photonic strength and the maximum possible group 
velocity |Vkij-'|. In Figurc[T]wc show the DOS in an empty 
fee crystal. As expected for a crystal with zero dielectric 
contrast, the calculated DOS is independent of dipole po- 
sition and orientation. In agreement with the well-known 
DOS in vacuum the calculated DOS increases paraboli- 
cally with frequency. Fluctuations around the parabola 
are due to binning noise associated with the finite k- 
space discretization. 

We have calculated the relative root-mean-square er- 
ror in the calculated density of states (DOS) for an 
fee 'empty' crystal averaged over all histogram bins 
in the frequency range < Lua/2TTc < 1 for several 
combinations of histogram binwidth and k-grid resolu- 
tion, as specified in the caption of Fig. [2] As predicted 
by Eq. (jlip . the deviation of the approximation from 
the analytic result is inversely proportional to the ratio 
Alo/ Ak. In most cases, one wants to calculate the LDOS 
with a given frequency resolution Auj, i.e. N(r, uj, ed)Aa;, 
to within a predetermined accuracy. For instance, calcu- 
lating the vacuum DOS for frequencies < uia/2-Kc < 1 
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with a desired absolute accuracy better than 0.01 • 
(1% of the vacuum DOS at Loa/iirc = 1) and a desired 
frequency resohition of Alj = 0.01(27rc/a), requires using 
Afc ^ Aw/O.Sc. The number of k-points corresponding 
to this wave vector samphng equals 2480 k-points in the 
irreducible wedge of the fee Brillouin zone, 59520 in the 
required half Brillouin zone, or cquivalcntly = 119040 
in the full Brillouin zone. Photonic crystals with nonzero 
index contrast cause a pronounced frequency structure 
of the LDOS. Due to the flattening of bands compared 
to the dispersion bands of vacuum-only crystals, the k- 
space integration itself is at least as accurate, assuming 
that the eigenfrequencies and the field-mode patterns are 
known with infinite accuracy. In practice, one needs to 
adjust the number of plane waves to obtain all eigen- 
frequencies to within the desired frequency resolution 
Aw. In our computations for fee crystals, we represented 
the k-space of a half of the first Brillouin zone by an 
equidistant grid consisting of Nk/2 ~ 145708 k-points. 
The frequency resolution of our LDOS histograms is 
Ao; = 0.01 • 27rc/a, and the spatial resolution of the 
LDOS is approximately a/40 judging from the maximum 
|G — G'l « 120/a involved in the expansion with Nq ~ 
725 plane waves in Eq. jl]). 



D. Computation time required for LDOS 

Calculating the LDOS in 3D periodic structures is a 
time-consuming task. The chosen degree of k-space dis- 
cretization {Nk/2 ~ 145708 k-points) and the number 
of plane-waves {Na = 725) are the result of a trade-off 
between the desired accuracy and tolerated duration of 
the calculations. Essentially the computation consists of 
solving for the lowest n eigenvalues and eigenvectors of 
a real symmetric 2Na x 2Nc matrix for each of the Nk 
k-points independently. To accomplish this calculation, 
we use the standard Matlab "eigs" implementation [43] 
of an implicitly restarted Arnoldi method (ARPACK) 
that takes approximately 2.2 seconds to find the low- 
est 20 eigenvalues and eigenvectors of a single Hermitian 
1500x1500 (i.e., Nq ~ 750 plane waves) matrix on a 
3 GHz Intel Pentium 4 processor, or on a 2.4 GHz In- 
tel Core Duo processor. For 145708 k-points the result- 
ing computation time is hence on the order of 90 hours 
(4 days). Since the Matlab ARPACK routine is already 
highly optimized, we do not expect that the computation 
time per k-point can be significantly reduced for LDOS 
calculations based on the standard H-field inverted ma- 
trix plane wave method. In terms of the number of plane 
waves Nc, the computation time scales as Nq. As in 
Ref. [44], the algorithm can be accelerated by realizing 
that the iterative eigensolver does not require the full 
matrix Ti. multiplying the vector u in Eq. [3 but rather 
a function that quickly computes the image 7iu of any 
trial vector u. Since calculating Hu repeatedly involves 
a slow matrix- vector multiplication with 77g_g', the al- 
gorithm is only accelerated for Nq > 1000. 




a/x 

Fig. 3. (Color onlinc)DOS per volume in an fee crys- 
tal consisting of spheres with e = 7.35 in a medium 
with e = 1.77 with a filling fraction of the spheres of 
25 vol %. The solid dotted curve represents calculations 
from Ref. [31]. Our result is plotted as a histogram (red). 



3. Comparison with previous results 

To test the computations, we compare our results with 
earlier reports. We have calculated the DOS and LDOS 
in an fee crystal consisting of dielectric spheres with ei 
— 7.35 (Ti02) in a medium with e2 = 1.77 (water) - 
the same structure as was analyzed in Refs [31,33]. The 
spheres occupy 25 vol% of the crystal. Figure [3] shows the 
total DOS in this photonic crystal calculated by us and 
by Busch and John [31] . We reproduce the earlier cal- 
culations of the total DOS: both results shown in Fig. [3] 
are in excellent agreement, with deviations less than 2% 
throughout the frequency range < Loajl-Kc < 1. In Fig- 
ure 2] (Media l)we demonstrate the LDOS in the same 
photonic crystal at a specific location: at a point equidis- 
tant from two nearest-neighbor spheres. In this calcu- 
lation, we used the same number of reciprocal-lattice 
vectors = 965 as in the only available benchmark 
paper by Wang et al. [33] that docs not contain sym- 
metry errors. We find that our calculations (empty cir- 
cles) are in good agreement up to a/A = 0.85 with the 
LDOS reported previously (solid curve): the deviations 
are smaller than 1%. Deviations at higher frequencies 
are either due to a difference in accuracy of k-space 
integration (frequency binning resolution and k-space 
sampling density), or to a difference in accuracy of the 
plane wave methods (eigenmodes and eigenfrequencies). 
Unfortunately, neither the accuracy nor the frequency 
resolution is specified in Ref. [33] . The k-space sampling 
density in our work is approximately twice the density 
specified in Ref. [33]. Based on the excellent agreement 
of our total DOS calculations with those of Busch and 
John [31], and on the fact that Wang ct al. used a k-space 
density and plane wave number comparable to that of 
Busch and John, it is unlikely that inaccuracy of the k- 
space integration in either calculation is the source of dis- 
crepancy. We therefore surmise that deviations are due 
to a difference in the method of evaluation of E„.k(r). 
Conversion of D to E by multiplication with l/e(r) in 
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Fig. 4. (Color online) Dipole-averaged LDOS in the same 
photonic crystal as in Fig. [3] at a position {j,j,0). Red 
histogram: our calculations (Media 1). Solid curve: re- 
sults from Ref. [33]. This relative LDOS is the ratio of 
the LDOS to that in vacuum at a/A = 0.495. 

real space as proposed in [33] can cause incorrect mode 
amplitudes due to Gibbs oscillations, as opposed to the 
Fourier space conversion using Eq. ([5]). In general our 
calculations confirm the result by Wang et. al. [33] that 
the LDOS is only correctly calculated by integration over 
half the Brillouin zone, rather than over the irreducible 
part as was incorrectly used in earlier reports. 

4. LDOS in Ti02 inverse opals 

In recent time-rcsolvcd experiments, enhanced and in- 
hibited emission rates were demonstrated for quantum 
dots embedded inside strongly photonic Ti02 inverse 
opals [19,21,35]. In the framework of these experiments, 
it is highly relevant to calculate the LDOS inside such 
inverse-opal photonic crystals, especially for the source 
positions occupied in experiments. We model the posi- 
tion dependence of the dielectric function e(r) as shown 
in Figure [5] This model assumes an infinite fee lattice 
of air spheres with radius r = 0.25'\/2a (a is the cubic 
lattice parameter). The spheres are covered by overlap- 
ping dielectric shells (e ~ 6.5) with outer radius 1.09r. 
Neighboring air spheres are connected by cylindrical win- 
dows of radius 0.4r. The resulting volume fraction of 
Ti02 is equal to about 10.7%. These structural param- 
eters are inferred from detailed characterization of the 
inverse opals using electron microscopy and small an- 
gle X-ray scattering [14, 15]. Moreover, the stop gaps in 
the photonic band structure (Fig. [6] (Media 2 and Me- 
dia 3)) calculated using this model agree well with re- 
flectivity measurements in the ranges of both the first- 
order (a/A = 0.7) and second-order Bragg diffraction 
(a/A = 1.2) [45]. 

Henceforth, we will consider the relative LDOS, which 
is the ratio of the LDOS in a photonic crystal to that in 
a homogeneous medium with the same volume-averaged 
dielectric function. This scaling is motivated by exper- 
imental practice, in which emission rate modifications 
are judged by normalizing measured rates to the emis- 
sion rate of the same emitter in crystals with much 
smaller lattice constant a [19,21,35]. In these reference 




Fig. 5. Rendering of the dielectric function in one fee unit 
cell that models the Ti02 inverse-opal structure in sec- 
tion[3| an fee lattice of air spheres of radius r = 0.25-\/2a 
with a being the cubic lattice parameter. The spheres 
are covered by shells with e = 6.5 and outer radius 
1.09r. Neighboring air spheres are connected by win- 
dows of radius 0.4r. The letters (a-d) indicate four dif- 
ferent positions at the Ti02-air: a = (1,0, 0)/(2v^), b = 
(1, 1, 2)/(4V3), c = (1, 1, 1)/(2V6) and d = (0.33,0.13,0) 
(points shown are symmetry-equivalents). The dash- 
dotted line shows the main diagonal of the cubic unit 
cell. 




Fig. 6. Left: Photonic band structure (Media 2)for the 
Ti02 inverse opal shown in Fig. [S) The grey rectangles 
indicate stopgaps in the PL direction and one stopgap 
in the PX direction for the inverse opal. The stopgaps 
result in the decreased DOS {right: circles, Media 3) at 
corresponding frequencies compared to the DOS in a ho- 
mogeneous medium with Uav = 1-27 {right: solid line). 
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Fig. 7. Relative LDOS in the inverse opal shown in Fig [5] 
at three different positions: (a, Media 4) r = (0,0,0) [the 
center of an air sphere, solid curve], r = i(f , 1, 1) [among 
three air spheres, dash-dotted curve] and r = (5, 0,0) 
[midway between two spheres along [1,0,0] direction, dot- 
ted curve]; (b. Media 5) r = i(l,l,0) [in the window 
between two spheres] projected on [1,1,0], [-1,1,0] and 
[0,0,1] directions shown by solid, dash-dotted and dot- 
ted curves, respectively. 

systems, emission frequencies correspond to the effective 
medium limit a/A ^ 0.5 quantified by an average index 
nav = >/eai; = 1-27 for the Ti02 crystals [19]. In units 
of 4/a^c, the LDOS in a homogeneous medium is equal 
to nav{a/X)^/3. In Figure [3i (Media 4) we plot the re- 
sulting LDOS at three positions in the unit cell: at r = 
(0,0,0), i(l, 1, 1) and (^,0,0). Due to the high symmetry 
of these points, the LDOS does not depend on the dipole 
orientation, as we verified explicitly. A first main obser- 
vation from Fig. [7^ (Media 4) is that the LDOS differs 
considerably between these three positions at all reduced 
frequencies. This observation illustrates the well-known 
strong dependence of the LDOS on position within the 
unit cell of photonic crystals [10,27,31]. 

A second main observation in Fig. [7^ (Media 4) is that 
the LDOS strongly varies with reduced frequency, reveal- 
ing troughs and peaks caused by the pseudogap near 
a/X = 0.7, which is related to P*-order stop gaps such 
as the L-gap. The effects of 2"''-order stopgaps appear 
beyond a/X > 1.15 [45]. In the middle of the air region 
at position r = (0,0,0) there is a sharp, factor-of-thrce 
enhancement at a/X w 1.35 within a narrow frequency 
range. This feature could be probed by resonant atoms 
infiltrated in the crystals [46]. At position r = (^,0,0) in 
an interstitial, the mode density has a broad trough near 
a/X = 1.25 that will lead to strongly inhibited emission. 

A third main observation is that at spatial positions 
with low symmetry, the LDOS clearly depends on the 
orientation of the transition dipole moment. Figure [TJd 
(Media 5) shows the frequency dependent LDOS for 
three perpendicular dipole orientations at a position in 
the center of a window that connects two neighboring 
air spheres (see Figure [5]). The LDOS differs for all three 
orientations, and is thus anisotropic. At low frequency 
{a/X = 0.3), the emission rate is highest for a dipole 
pointing in the (1,1,0) direction, with increasing fre- 
quency the highest rate shifts to the (0, 0, 1) and then to 
the (—1,1,0) orientation. We emphasize that for emit- 
ters with fixed or slowly varying dipole orientations, such 




position [ria] position [ria] 



Fig. 8. Relative LDOS at four key frequencies, a/X = 
0.535, 0.725, 0.865 and 1.295, in the inverse opal as a 
function of position r from (0,0,0) in the [1,1,1] direction. 
The hatched boxes indicate the position of the dielectric 
Ti02 shell. The LDOS is projected on two dipole ori- 
entations: (a) [1,1,1] perpendicular to the dielectric-air 
interface, and (b) [-1,1,0] parallel to the interface. The 
LDOS projected on the [-1,-1,2] and [-1,1,0] directions is 
equal. For r/ae[^,VS\ the LDOS is mirror-symmetric 
to that in the region [0, -^]. 



as dye molecules or quantum dots on solid interfaces, the 
emission rate is determined by the optical modes that are 
projected on the dipole orientation. Therefore, knowl- 
edge of the projected LDOS is important for controlling 
spontaneous-emission rates as well as for interpreting the 
data from experiments on emitters in photonic metama- 
terials. 

A fourth main observation from Fig. [7^ and b (Media 
4 and 5) is that at low frequencies a/X < 0.5, the rela- 
tive LDOS hardly varies with frequency, which means 
that the mode density is proportional to uj'^, as in homo- 
geneous media. Interestingly, there is a clear dependence 
of the LDOS on both the position and the dipole orienta- 
tion even at these low frequencies, i.e., long wavelengths 
relative to the crystal periodicity. The reason for these 
effects is that the photonic Bloch modes exhibit local 
variations of the electric field related to local variations 
of the dielectric function in order to satisfy the conti- 
nuity equations at dielectric boundaries for the parallel 
E or perpendicular D field, respectively [34,47]. Conse- 
quently, the LDOS strongly varies on length scales much 
less than the wavelength, i.e., even in electrostatic or ef- 
fective medium limit. While such behavior may appear 
surprising, its origin in electrostatic depolarization ef- 
fects has been discussed before [47, 48] . 

To gain more insight in the spatial dependence of the 
LDOS in the inverse opals, we have performed calcula- 
tions for dipolcs positioned along a characteristic axis in 
the unit cell. Figure [8] shows the LDOS at four key fre- 
quencies for dipoles placed on the body diagonal of the 
cubic unit cell, that is, from (0, 0, 0) in the [1, 1, 1] direc- 
tion. The LDOS has clear maxima and minima along the 
diagonal, varying by more than 10 x. Figure [51^ a) shows 
that for a dipole oriented in the [1,1,1] direction (per- 
pendicular to the dielectric-air interface), the LDOS is 
strongly (> 5x) suppressed near the dielectric Ti02 shell 
{r/a ~ 0.353) at all frequencies. In contrast, no strong 
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Fig. 9. Relative LDOS in the inverse opal at four different 
positions on the Ti02-air interface shown in Fig. [5j At 
each position the LDOS is projected on three mutually 
orthogonal dipole orientations e^. (a, Media 6) point a 
for Gd = [1,0,0] and [0,1,0]. LDOS at = [0,0,1] and 
[0,1,0] is the same, (b. Media 7) point b for e^; = [1,1,2], 
[-1,1,0] and [-1,-1,1]. (c. Media 8) point c for e<j = [1,1,1] 
and [-1,1,0]. LDOS at = [-1,-1,2] is equal to that at 
ed = [-1,1,0]. (d. Media 9) point d for = [0.33,0.13,0], 
[-0.13,0,0.33] and [0,0,1]. 



suppression or enhancement occurs for dipole orienta- 
tions parallel to the dielectric interface, see Fig. [H^b). 
The strongly differing LDOS for dipoles near the dielec- 
tric rationalize the broad distributions of emission rates 
that were recently observed for quantum dots in inverse 
opals [35] , see below. Enhancements and inhibitions also 
occur in the air regions: the LDOS is enhanced at r/a ~ 
0.28 and 0.57 by up to 2.5 to 3 times, respectively, see 
FigureO^a). At point r/a = 1/^/3 « 0.57 that lies in the 
(111) lattice plane in the air region, the LDOS is inhib- 
ited at all frequencies for the dipole orientations [-1,1,0] 
and [-1,-1,2] that are perpendicular to the (111) plane 
(see Fig. IHb)). Finally, for dipoles parallel to the body 
diagonal (Figure [8]Ja)) the mode density shows pseudo- 
oscillatory behavior on length scales much less than the 
wavelength, e.g., a period of 0.2a at frequency 0.535a/A 
(period corresponds to 1/10*'' of a wavelength). This ob- 
servation confirms that photonic crystals are bona fide 
metamaterials where optical properties strongly vary on 
length scales much less than the wavelength. 

In recent time-resolved experiments [21,35], emission 
from quantum dot light sources distributed at the in- 
ternal Ti02-air interfaces of the inverse opals was in- 
vestigated. To analyze the experimental data, we have 
calculated the LDOS at several symmetry-inequivalent 
positions on the Ti02 shells: at points a, b, c and d 
(see Fig. [5]) for three mutually orthogonal orientations 
of the emitting dipole, where the first orientation is cho- 
sen along the vector pointing from (0,0,0) toward the 
corresponding point. Figures [9]Ja) through [9]Jd) (Media 
6 through 9) show that at all these positions the LDOS 



strongly varies with reduced frequency and position, as 
expected, and also with orientation of the dipole. In 
broad terms, for a dipole parallel to the interface the 
LDOS is near 1 at low frequency, increases to a peak 
enhancement of 2.5 x at a/X ~ 1.22 before strongly 
varying at high frequencies. For reference, the frequency 
a/A = 1.22 is in the range where a band gap is expected 
for more strongly interacting crystals. The plots also re- 
veal that for dipoles perpendicular to the Ti02-air in- 
terface, the LDOS is strongly inhibited (more than lOx) 
over broad frequency ranges. For instance. Figure [DJa) 
shows that the emission rate for a dipole perpendicular 
to the interface is 16-fold inhibited to a level of 0.06, with 
a maximum of 0.22 (5-fold inhibition) at a/A = 1.22. 
It is striking that the strong inhibition occurs over a 
broad frequency range from 0.3 to 1.6, i.e., more than 
two octaves in frequency. While the inhibition is not 
complete, the bandwidth is much larger than the maxi- 
mum bandwidth of 12% for a full 3D bandgap in inverse 
opals [49] and far exceeds the bandwidth of the 2D gap 
for TE modes in membrane photonic crystals, that allows 
a seven-fold inhibition [26, 30] over a 30% bandwidth. 

In the frequency range up to a/X = 1.1, the depen- 
dence of the LDOS on frequency is quite similar at all 
the positions and dipole orientations. This remarkable 
result agrees with the experimental observation from 
Ref. [35]: the complex decay curves of light sources in 
the inverse opals are described by one and the same func- 
tional shape (log-normal) of the decay-rate distribution 
for all reduced frequencies studied. 

For the interpretation of the time-resolved experi- 
ments on ensembles of emitters in the inverse-opal pho- 
tonic crystals, the results presented above mean that 

1. the decay rate of an individual emitter is determined 
by its frequency, position and also by the orientation 
of its transition dipole in the photonic crystal; 

2. the measured spontaneous-emission decay depends 
on how the emitters are distributed in the crystal; 

3. even in the low-frequency regime, ensemble 
measurements will reveal non-exponentional decay 
curves; 

4. the similar shape of the reduced-frequency de- 
pendence of the LDOS allows modeling of the 
non-exponentional decay curves with a single type 
of decay-rate distributions. In other words, one 
distribution function can be successfully used to 
model the multi-exponentional decay curves meas- 
ured from crystals with different lattice parameters. 

5. LDOS in silicon inverse opals 

A complete inhibtion of spontaneous emission may only 
be achieved in photonic crystals with a 3D photonic 
band gap. Therefore, there has been much effort to fabri- 
cate inverse opals from silicon rather than titania, since 
the higher index of silicon allows for a photonic band 
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Fig. 10. Left: Photonic band structure (Media 10) for 
an inverse opal from silicon (e = 11.9). Right: The to- 
tal DOS in the Si inverse opal (Media 11). The DOS is 
strongly depleted for frequencies near FL and FX stop- 
gaps (grey rectangles). A photonic band gap (grey bar) 
occurs between bands 8 and 9, as also reflected in the 
vanishing DOS. 
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Fig. 11. Relative LDOS in a Si inverse opal at: (a, Me- 
dia 12) r = (0,0,0) [the center of an air sphere, solid 
curve], r = ^(1,1,1) [among three air spheres, dash- 
dotted curve] and r = (^,0,0) [midway between two 
spheres along [1,0,0] direction, dotted curve]; (b. Media 
13)r=-|(l,l,0)[in the window between two air spheres] 
projected on [1,1,0], [-1,1,0] and [0,0,1] directions shown 
by solid, dash-dotted and dotted curves, respectively. 



gap [31,37]. For the calculation of the LDOS in such Si 
inverse opals, the dielectric function e(r) was modeled 
similarly to that in the inverse opals shown in Fig. O 
From SEM observations we inferred the following struc- 
tural parameters [17, 18]: the outer radius of the over- 
lapping dielectric shells with e = 11.9 is about 1.15r 
(where r = 0.25\/2 a is the air-sphere radius, a is the 
lattice parameter). The cylindrical windows connecting 
neighboring air spheres have a radius of 0.2r. The larger 
outer radius and the smaller window size compared to 
the Ti02 inverse opals arc commensurate with a higher 
volume fraction of about 23% Si [18]. 

The band structure and total DOS for this system are 
shown in Fig. [TO] (Media 10 and Media 11). The DOS is 
strongly depleted in a pseudo-gap between the 2"'' and 
3'"'' bands [36], at frequencies near 0.55. Near frequency 
0.85, both the band structures and the DOS reveal a 
3D photonic band gap of relative width Acj/tu « 3 % 



between the 8*'* and 9*'* bands [31]. Compared to the 
Ti02 structure, both the lowest-order L-gap and the 8*'' 
and 9*'* bands are shifted to lower reduced frequencies. 
This shift is due to the higher effective refractive index of 
the Si inverse opals (riav = 1.88), on account of a higher 
index of the backbone and a higher filling fraction. 

Figure \TT\a) (Media 12) presents the relative LDOS 
at three high-symmetry positions in the unit cell r = 
(0,0,0), 1(1,1,1) and (5,0,0). The LDOS varies much 
more strongly with frequency than in Ti02 inverse opals, 
as a result of the larger dielectric contrast that leads to 
strongly modified dispersion relations and Bloch mode 
profiles. While the maxima up to 3.2 are not much higher 
than in Ti02 inverse opals, the minima in the mode den- 
sity arc reduced and the slopes arc steeper. As expected, 
the LDOS is completely inhibited at all positions in the 
frequency range of the band gap. Previously, it has been 
suggested that the LDOS could be inhibited at salient 
positions in the unit cell for frequencies outside a com- 
plete gap. While Figure [TlT a) (Media 12) reveals that 
the mode density is strongly reduced above the band 
gap, e.g., at r = (i,0,0), it is not truly inhibited. In 
fact, in the course of our study, we have not encoun- 
tered any "sweet spots" where the LDOS is completely 
inhibited at frequencies outside a 3D photonic band gap. 

Figure [TT^b) (Media 13) shows the LDOS at a lower 
symmetry position r = -1(1, 1, 0), in the window between 
two nearest-neighbor air spheres. The mode density has 
been calculated for three perpendicular dipole orienta- 
tions. Figure [lUb) (Media 13) shows that the LDOS is 
highly anisotropic since it differs for all three orienta- 
tions: at low frequencies (a/A < 0.5), the mode density 
is highest for the (—1, 1,0) orientation, intermediate for 
the (0,0,1) orientation, and lowest for the (1,1,0) ori- 
entation. With increasing frequency, the highest LDOS 
also changes to other orientations, with the (0, 0, 1) ori- 
entation having the highest LDOS above the pseudo-gap 
and even a narrow peak at frequency 0.95. Therefore, if 
it is possible to orient a quantum emitter with its dipole 
parallel to (0,0, 1), this frequency range is conducive to 
strong emission enhancement and perhaps even QED ef- 
fects beyond weak-coupling. Interestingly, these frequen- 
cies are slightly higher than the upper edge of the band 
gap where strong coupling effects have recently been dis- 
cussed [12]. 

6. Conclusions 

We have performed intensive calculations of the local 
density of states in Ti02 and Si inverse opals with exper- 
imentally relevant structural parameters. Since conflict- 
ing and incorrect reports have appeared on the LDOS in 
photonic crystals, we have set out to validate our method 
of choice, i.e., the H-field plane- wave expansion method. 
This validation relied on comparison to literature results, 
on the explicit verification of required symmetries that 
previous reports failed to satisfy, and on quantitative 
considerations of resolution and accuracy. Results for 
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each structure are made available for other workers in 
the field, both as benchmarks and for comparison with 
experimental data. With the help of these computations 
we have obtained quantitative insight in the LDOS rel- 
evant for time-resolved ensemble fluorescence measure- 
ments on photonic crystals, such as obtained in recent 
experimental work [35] . The results of our numerical cal- 
culations reveal a surprisingly strong dependence of the 
LDOS on the orientation of the emitting dipoles. 
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